Nonlinear screening and stopping power in two-dimensional electron gases 



E. Zaremba 

Department of Physics, Queen's University, Kingston, Ontario, 
Canada K7L 3N6 and Donostia International Physics Center (DIPC), 
Paseo Manuel de Lardizabal, no. 4> 20018 San Sebatian, Basque Country, Spain 

I. Nagy 

Department of Theoretical Physics, Institute of Physics, 
Technical University of Budapest, H-1521 Budapest, Hungary, 
and Donostia International Physics Center (DIPC), Paseo Manuel de Lardizabal, 
no. 4- 20018 San Sebatian, Basque Country, Spain 

P. M. Echenique 
Departamento de Fisica de Materiales, Facultad de Quimica, 
Universidad del Pais Vasco and Centro Mixta CSIC-UPV/EHU, 
Apdo. 1012, 20018 San Sebastian, Basque Country, Spain 
(Dated: February 2, 2008) 

We have used density functional theory to study the nonlinear screening properties of a two- 
dimensional (2D) electron gas. In particular, we consider the screening of an external static point 
charge of magnitude Z as a function of the distance of the charge from the plane of the gas. The self- 
consistent screening potentials are then used to determine the 2D stopping power in the low velocity 
limit based on the momentum transfer cross-section. Calculations as a function of Z establish the 
limits of validity of linear and quadratic response theory calculations, and show that nonlinear 
screening theory already provides significant corrections in the case of protons. In contrast to the 
3D situation, we find that the nonlinearly screened potential supports a bound state even in the high 
density limit. This behaviour is elucidated with the derivation of a high density screening theorem 
which proves that the screening charge can be calculated perturbatively in the high density limit for 
arbitrary dimensions. However, the theorem has particularly interesting implications in 2D where, 
contrary to expectations, we find that perturbation theory remains valid even when the perturbing 
potential supports bound states. 

PACS numbers: PACS Numbers: 73.20.-r, 73.20.Hb 



I. INTRODUCTION 

The screening response is a fundamental property of an electron gas in arbitrary dimensions. The situation in two 
dimensions is of particular interest because of the possible realization of quasi-two-dimensional systems in semiconduc- 
tor heterostructuresi, image and band-gap surface states at metal surfaces^, electrons at the surface of liquid helium*, 
metallic overlayers on insulating substrates and layered materials^. Examples of problems in which an understanding 
of the screening of an external charge is important include impurity-limited electron transport in two-dimensional 
electron gasesi, the electronic structure of intrinsic or photo-induced defects in semiconductor heterostructures^ and 
dynamic interactions with moving charges^. In each of these cases, the strong interaction of mobile electrons with 
the charged impurity leads to significant modifications of the local electronic structure. 

One of the first attempts to deal with 2D screening was the work of Stern and Howard^ which was motivated by 
measurements of electron mobilities in Si inversion layers. To perform explicit calculations they considered a model 
in which the two-dimensional electron gas (2DEG) is represented as a sheet of zero thickness, the ideal 2D limit, 
with a charged impurity situated a distance d from the sheet. The screened impurity potential was then obtained 
within a Thomas-Fermi approximation and was used to investigate the possibility of impurity bound states. However 
no attempt was made to perform fully self-consistent calculations of the screening charge density and the screened 
potential. This was taken up later in a series of papers by Vinter— for a Si inversion layer. This work is especially 
notable for its attempt to realistically account for the underlying bandstructure of the semiconductor and the finite 
thickness of the 2DEG. A detailed discussion of screening and bound state formation was presented for this particular 
heterostructure and the results were shown to be consistent with observed impurity-limited mobilities. 

Also of interest is the dependence of screening on the dimensionality of the system. That important differences 
might arise is indicated by the well known fact that any purely attractive potential in 2D always has at least one 
bound stateAP, in sharp contrast with the situation in 3D. Extending previous nonlinear screening calculations from 
3Dii to 2D and exploring these possible differences is part of the motivation for the work described in this paper. To 
do this we perform self-consistent nonlinear screening calculations for an ideal 2DEG within the context of density 
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functional theory. Our studies as a function of density reveal some interesting properties of the screening in 2D that 
had not been appreciated before. The interpretation of these numerical results is facilitated by the proof of a general 
theorem regarding the nature of screening in the high density limit. A preliminary account of some of this work was 
presented earlier—. 

As a specific application of the results we consider the low velocity stopping power for a heavy projectile moving 
parallel to the plane of the 2DEG. This problem has been treated previously in both linear and nonlinear response 
approximations. Of the linear theories, we mention the early work of Fetter^ using a hydrodynamic model which was 
followed by a calculation based on the random phase approximation (RPA) by Horing et aim*. In this latter work a 
finite distance, d, between the projectile and the plane of the gas was considered. The case of a projectile moving 
in the plane of the gas has also been studied, not only within the RPAi£, but also with the inclusion of local field 
corrections- 16 and for a gas at finite temperatures^. Of course the interaction cannot always be considered as weak 
and a fully nonlinear theory of the screening is in general necessary. One step in this direction is the RPA quadratic 
response formulation 18 which was used by Bergara et aZ— in 2D to provide the lowest order nonlinear correction to 
the stopping power. Although this theory is perturbative and therefore limited in its range of validity, it has the merit 
that the nonlinear correction is determined for arbitrary projectile velocities. 

The stopping power at low projectile velocities can also be formulated in terms of the scattering of electrons from the 
screened potential. For qualitative purposes it is sometimes useful to consider model potentials 2 ^ whose arbitrariness 
can be limited by imposing physical constraints such as the Friedel sum rule. This was the approach taken by Wang 
and Ma^i to determine proton and antiproton stopping powers in a 2DEG. A somewhat more fundamental approach 
based on the nonlinear screening theory of Sjolander and Stotli 2 ^ was used by Krakovsky and Percus 2 ^, for negatively 
charged projectiles. Their results for this case are consistent with the model potential results of Wang and Ma2i. 
None of these calculations, however, are truly self-consistent. Our results based on self-consistent density functional 
theory calculations eliminate this source of uncertainty in the determination of 2D stopping powers. 



II. SELF-CONSISTENT SCREENING 



The electrons of our ideal 2DEG are assumed to have an isotropic effective mass m* and move in the presence of 
a uniform neutralizing positive background. We also imagine the entire 2D system to be immersed in an extended 
dielectric with permittivity e as might occur in a semiconductor heterostructure. Other dielectric arrangements can 
also be envisaged, for example the surface of liquid helium, but these will not be considered here. All lengths will 
be expressed in units of the effective Bohr radius ao = eh 2 / m* e 2 , and all energies in units of the effective Hartree, 
H = e 2 /eclq. The mean density of the gas, rto, is characterized by the density parameter r s = 1/ \/im$. 

Our objective is to determine the nonlinear screening of a stationary point charge, Z, located a distance d from the 
plane of the 2D gas. The methodology follows that used in earlier calculations for a three-dimensional gasiii 2 ^ and is 
suitable for both positive and negative charges. This screening response is determined by solving self-consistently the 
two-dimensional Kohn-Sham equations 

- iv 2 V^(r) + A«eff(r)^(r) = E^(r) (1) 

where the effective potential is given by 

Auoff (r) = v ext {r) + Av H (r) + Av xc {r) . (2) 
The external potential v ext (r) in the plane of the gas is — Z/^/d 2 + r 2 , and Avh(y) is the Hartree potential 

AMD = / *>f$ (3) 

due to the electronic screening density, Ara(r) = n(r) — uq. The change in the exchange-correlation potential, 
Av xc (r) — v xc [n(r)] — v xc [no], is defined in the local density approximation (LDA) using the parameterization of the 
2D exchange-correlation energy given by Tanatar and Ceperley— 
The total screening density is given by 

An(r) = \Mr)\ 2 + £ [IV>,(r)| 2 - |V 4 °(r)| 2 ] (4) 

b i 

where the first sum extends over all bound states of the effective potential, and the second extends over all occupied 
continuum states up to the Fermi level Ep. We assume that each spatial orbital is doubly-occupied for spin. The 
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free-particle solutions ^>°(r) are obtained in the absense of the external potential. Due to the cylindrical symmetry 
of the problem, the solutions of Eq. have the form ipi(r) = e lm ^R nm {r) where the angular momentum quantum 
number m takes on integral values and n distinguishes the different radial solution o?i 2 9 . The continuum states behave 
asymptotically as Rkm(f) ~ r" 1 / 2 cos(fcr — ^\m\n — in + r) m ) where k — \/2E is the wavevector and rj m {k) is the 2D 
scattering phase shift. For the free-particle solutions, r\ m — 0. These scattering phase shifts are related to the total 
screening density according to the 2D Friedel sum rule (FSR)£ 

2 °° 

Zfsr = ~ T] m (E F ). (5) 

m=— oo 

At self-consistency, this quantity must equal the charge Z of the external impurity. 

The solutions of Eq. (JTJ are obtained using standard numerical techniques. The calculation of Aujj(r), the only 
nontrivial step in determining the effective potential, is performed by making use of a fast fourier transform technique 
to go successively between real and wavevector space. Details of the method are given in an Appendix. 

III. NUMERICAL RESULTS 

Self-consistent calculations were performed for a range of densities and impurity charges Z , including negative values. 
In addition, calculations were done with and without exchange-correlation effects included in order to compare with 
earlier linear calculations which were performed within the random phase approximation (RPA). 

A. Z = 1, d = Q 

To begin we show results for a charge Z = 1 in the plane of the 2DEG (d = 0). Fig. 1 presents the self-consistent 
effective potential A.v a g for r s values ranging between 1 and 10 in steps of 1. This range spans low to high densities 
and in particular includes densities of experimental interest. Also shown for comparison is the Thomas-Fermi (TF) 
screened potential which, as will be explained shortly is the r s — > limit of the DFT calculations. As such, we see 
that the potential for r s = 1 is already quite close to the TF limit. 

The variation with r s is systematic, evolving with increasing r s into a potential with a strongly repulsive region 
beyond r ~ 0.5 a.u. This behaviour can be explained in terms of the underlying electronic structure of the screened 
impurity. For all the densities shown, the screened potential supports one bound state which is occupied by two 
electrons. Since the total screening density must integrate to unity according to the FSR in Eq. JSJ, the continuum 
states must themselves contribute a total charge of +1 in order to compensate for the overscreening provided by the 
two bound electrons. This is illustrated in Fig. 2 which shows the bound and continuum densities for two representative 
r s values, r s = 2 and 10. It can be seen that the bound state distribution is very similar despite the very different 
background densities. This is consistent with the similarity of the inner region of the potentials in Fig. 1 which is 
primarily responsible for the shape of the bound state. At large distances from the impurity the densities in Fig. 2 
exhibit characteristic Friedel oscillations with wavelength ir/kF- 

The negative portion of the screening densities in Fig. 2 corresponds to a local charge density which is positive 
and provides the necessary compensation of the negative bound state electronic charge. However, the distribution in 
space of this compensating charge is quite different in the two cases and leads to the different behaviour of the r s = 2 
and r s = 10 potentials for r > 0.5 a.u. The form of the r s = 10 continuum density is particularly interesting. In the 
range < r ~ 10 a.u., An cont (r) has approximately a constant value of —no (hence rAn cont (r) behaves linearly as 
seen in Fig. 2). In other words, the Z = 1 impurity with 2 bound electrons can be viewed as an H ion which sits at 
the center of a positively charged disc of radius R ~ r s . Returning to Fig. 1, we can now understand the repulsive 
part of the potentials seen for large r s . The LDA potential of a 2D H~ ion itself has a +l/r Coulomb tail arising 
from the Hartree potential. In the 2DEG environment, this H - potential is screened by the positive background and 
as a result, the potential goes to zero at approximately R ~ r s . 

Another interesting quantity is the bound state energy shown in Fig. 3 as a function of r s . The solid curve is the 
result of the full DFT calculation including xc, while the dashed curve is the corresponding result obtained when Av xc 
is set to zero. We refer to the latter as the Hartree approximation. The Hartree result is monotonically decreasing 
with decreasing r s and goes to a limiting value of Eq = —0.2862 H at r s = which, as we shall see, is the binding 
energy for the TF screened potential. The xc result likewise decreases with decreasing r s , but reaches a mininum near 
r s ~ 0.5 and then increases to the same limiting value at r s — as in the Hartree approximation. The more negative 
xc eigenvalue reflects the stronger binding due to the xc potential which is attractive in the vicinity of the positive 
impurity as a result of the pile-up of the electron screening density. To illustrate this we compare the self-consistent 
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potentials with and without xc for r s = 2 in Fig. 4. It can be seen from this figure that Av e g with xc is more attractive 
in the core region which is important for the determination of the bound state. Associated with this is a relative 
phase shift of the long range Friedel oscillations between the two cases. 

The decreasing trend in the xc eigenvalue is consistent with the behaviour of the potentials in Fig. 1 but is opposite 
to what is found in the analogous 3D calculations— There the bound state energy increases with decreasing r s , that is, 
the bound state becomes shallower with increasing density. The explanation for the 3D behaviour is that the screening 
of the impurity potential is more effective with increasing density, and as a result, the bound state eventually ceases 
to exist beyond a certain critical density. This behaviour is consistent with the expected applicability of perturbation 
theory in the high density limit when the Fermi energy Ep is much larger than the magnitude of the screened potential 
| Aw c ff (r)J— . If one were to use the same reasoning in 2D, the screening density in wavevector space would be given 
by the linear response result (Sec. IV provides further details) 

M?) = Xo(q)&v cS (q) , (6) 

where Xo(l) is the Fourier transform of the static noninteracting density response function of the 2D gas. Within 
linear response theory, the screening density arises from the perturbation of the plane wave states of the uniform gas 
and there cannot be a bound state contribution. However, as we have just seen, bound states do persist in 2D even 
in the high density limit. This is a peculiarity of 2D and is associated with the fact that a purely attractive potential 
always has at least one bound state regardless of the strength of the interaction^. This effect is obviously missed 
by linear response theory and would seem to invalidate a perturbative approach. Nevertheless, it turns out that the 
results of linear response theory are indeed correct even though the theory neglects bound states entirely. We defer 
discussion of this point to the following Section. 

The opposite limit r s — > oo is also of interest. In the Hartree approximation, the bound state eigenvalue approaches 
zero near r s — 8. Whether a very shallow bound state persists beyond this value could not be confirmed due to 
difficulties in obtaining converged self-consistent solutions. This was less of a problem with the inclusion of exchange 
and correlation since these additional interactions stabilize the bound state relative to the Hartree calculation. In 
the range 8 < r s < 12 the bound state eigenvalue could be fit approximately to the expression Eq ~ 0.023 — 0.7/r s 
H. Extrapolation of the numerical data beyond r s — 12 suggests that the bound state eigenvalue goes to zero at 
some large, but finite, r s value. However, it was not possible to confirm this since it was difficult to obtain converged 
solutions for large values of r s . Thus we cannot state unequivocally what the large r s behaviour is, and in particular, 
whether or not a free H~ ion is stable in 2D within the LDA. In 3D it is knownSS that the H~ ion is not stable within 
the LDA and our data suggest that this might also be the case in 2D. This remains an interesting question to explore 
in the future. 

As a final comment about the eigenvalue obtained with xc, we note that the H ion is stabilized in the range 
studied by the fact that it sits at the center of a positively charged disc of radius R~ r s . The electrostatic potential 
due to this charge has the value —2/r s H at the centre of the disc relative to infinity. The depletion of charge 
also gives rise to a shift in the xc potential of — v xc (no) which is positive and therefore a destabilizing effect. For 
large r s , —v xc (no) ~ +1.6/r s H which, when combined with the electrostatic shift, gives a potential shift of about 
— 0A/r s H . It is clear that this represents a significant contribution to the xc eigenvalue at the larger r s values 
we have considered. However the existance of a bound state cannot be determined without a detailed knowledge 
of the self-consistent potential near the impurity. In the case of the Hartree approximation, the stabilizing effect 
of the electrostatic potential of the positively charged disc (— 2/r s H) is clearly insufficient to compensate for the 
destabilizing effect of the Hartree interaction between the two bound electrons in the large r s limit. 

B. Z = -1, d = 

We consider next the situation of a negatively charged impurity such as an antiproton or acceptor state in a 
semiconductor—. Fig. 5 illustrates the variation of the screening density with ?' s . The results found here in 2D 
are similar to those found previously in 3D 2 4. The negative impurity repels electrons from its vicinity and roughly 
speaking exposes a positively charged disc of radius R ~ r s which neutralizes the impurity. In other words, the 
impurity sits at the center of a hole in the electron gas. The similarity to the situation for a positive impurity is 
emphasized in Fig. 6 where we superpose the screening densities for Z = 1 and Z = — 1 at r s = 10. The Z = — 1 
screening density is in fact very similar to the continuum part of the Z = I screening cloud. Thus, at sufficiently 
large r s , the H ion for Z = 1 effectively acts as an external Z = — 1 impurity as far as the rest of the electron gas is 
concerned. This emphasizes that the H~ configuration exists as a well-defined entity in this range of r s values. 

In all cases we find that the screened Z = — 1 potential does not support any bound states. We mention this since it 
has previously been claimed^ that the introduction of a negative test charge into a 2D gas could give rise to potentials 
which would bind an electron (or rather, a second negative test charge). To make contact with this earlier work we 
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compare in Fig. 7 the nonlinearly screened potentials with those obtained on the basis of linear response theory. One 
curve shows the screened potential as obtained in the RPA. This potential is defined in Fourier space as 

V RPA () = Kg) = Wg (7) 
V sc W e RPA^ e RPA( q ) > V) 

where the RPA dielectric function is given by 31 

e RPA (q) = l + v(q) Xo ( q ). (8) 

At this level of approximation, only Hartree interactions are relevant and the screened potential defined above is 
the potential experienced by a Z = —1 test charge moving in the vicinity of the impurity. A second curve shows 
the linearly screened potential, v^ c FC = v(q)/e LFC , obtained when local field corrections (LFC) are included in the 
definition of the dielectric function. Within the LDA, this is given by 

1 =1 %)Xo(g) (g) 



e LFC (q) l + (v(q) + v' xc ) X o(q) 

where v' xc (no) — dv xc (no) / dno is the LDA local field correction, i.e., in standard notation, v(q)G(q) — —v' xc (no)- The 
potential v^f c is simply the Hartree potential obtained with a screening charge density that includes the effects of 
LFC. As can be seen in Fig. 7 for r s — 4, this potential has a large negative region in real space and supports a bound 
state for a unit negative test charge of one electron mass. The authors of Ref. 3 — then argue that the potential also 
binds an electron in the 2DEG and that this binding may be a relevant pairing mechanism. However this conclusion is 
invalid for several reasons. First, the Hartree potential for the nonlinearly screened impurity (the short-dashed curve 
in Fig. 7) has a much shallower attractive portion as compared to the corresponding result of the linear calculation 
(long-dashed curve). A similar result was found analytically for a point charged screened by a uniformly charged 
disk22. But more importantly, an electron, as opposed to a negative test charge, also feels the effect of the induced xc 
potential. With this contribution included, the full self-consistent nonlinear potential Av c s defined in Eq. has a 
very shallow attractive part, and there is no tendency for bound state formation, as confirmed numerically. At a more 
fundamental level, the question of pairing involves the effective interaction between pairs of electrons at the Fermi 
level which typically have a high relative velocity. From this point of view it is unclear to what extent a statically 
screened potential for a test charge can be used as an estimate of the pairing interaction. 



C. Z = l, d^O 



We next consider the effect of moving the impurity charge out of the plane of the 2DEG. This is the situation 
corresponding to an external charge incident on the 2DEG from the outside, or a remote ionized impurity in a 
semiconductor heterostructurei^ 

Fig. 8 shows the screening density for a few values of d. It can be seen that the amplitude of the induced density 
decreases rapidly with increasing d and the density becomes more spread out along the plane. The reason for this 
behaviour is that the potential loses its Coulomb singularity as soon as the impurity is out of the plane, and therefore 
becomes smoother and weaker. In Fig. 9 we show the nonlinear screening density for the case r s = 2 and d = 2 
a.u. together with the bound state and continuum contributions. The bound state in this case is very shallow and 
extended, in fact more extended than the total density itself. This is possible since the continuum contribution to the 
density has a similar extent but is of opposite sign, so that there is a cancellation of the densities in the asymptotic 
region. 

Also shown in Fig. 9 is the screening density as determined by linear response theory. Despite the fact that the 
nonlinear screening response has a bound state contribution, the linear response result is seen to agree very well. The 
reason for this agreement is not obvious and will be explained in detail in Sec. IV. In addition, the density is already 
fairly close to the classical image theory result (d/2n)(r 2 + d 2 )~ 3 / 2 . This expression for d — 2a$ gives a density at the 
origin of 0.04 a.u. as compared to the nonlinear result of 0.035 a.u. With further increases in d the nonlinear density 
rapidly approaches the classical screening charge density. 



D. Stopping Power 



As discussed in the Introduction, there have been a number of calculations of 2D stopping powers. Our purpose 
here is to present results which are based on the full nonlinear screening charge densities and potentials. Within the 
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so-called kinetic theory framework, the stopping power at low velocities for a projectile moving with velocity v in the 
plane of the 2DEG is given by the expression 3 ^ 

S = n Q vv F a tr (E F ) , (10) 
where at r {E F ) is the momentum-transfer cross section defined in terms of the scattering phase shifts 

4 00 

<Jtr{E F ) = sin - 2 [nm(E F ) - Vm+l(E F )} . (11) 

VF m=0 

To leading order in the velocity, it is sufficient to determine the scattering phase shifts using the static nonlinearly 
screened potentials calculated in the present paper. 

In Fig. 10 we show the stopping power as a function of the projectile charge Z for d = and r s — 2. For small Z, 
S has the expansion S = S\Z 2 + S2Z 3 + . . . where the first two terms are the linear and quadratic response results, 
respectively. To emphasize the deviations from linear response, we present the results in the form S/(vZ 2 ). In this 
representation, the stopping power including the quadratic response correction appears as a straight line with slope 
S2/V. This correction was previously calculated within quadratic RPA— and is shown in Fig. 10 as the straight line. 
The fact that the line is tangent to the nonlinear Hartree curve at Z = shows that the dynamic quadratic response 
formulation is consistent with the kinetic theory approach for small values of Z and for low velocities. However, the 
validity of the quadratic response result is mainly limited to negative charges down to Z ~ — 1 and breaks down 
completely for Z ~ 1. It is thus apparent that nonlinearities are generally very important. Furthermore, it can be 
seen that the inclusion of xc leads to a large enhancement of the stopping power in the range — 1 < Z < 1, even to a 
greater extent than found in 3D. 9* This emphasizes that properties associated with the scattering of continuum states 
will be strongly affected by the xc part of the self-consistent potential. 

The behaviour of the stopping power as a function of density is also of interest and in Fig. 11 we show r s S/v vs. r s 
for Z = 1. The quantity S/v has the physical interpretation of a friction coefficients^ 3 , and is actually a monotonically 
decreasing function of r s with a finite limiting value of ttZ 2 for r s — > 0. This result follows^ from the 2D momentum- 
transfer cross section ot r {E F ) = Z (2ir / v F ) tanh(ir Z / v F ) which is obtained when scattering from a bare Coulomb 
potential, —Z/r, is considered. This is the exact result in the high density limit and simply reflects the fact that 
screening is not important for high energy scattering. As a result of this finite limiting value, the quantity r s S plotted 
in Fig. 11 shows a maximum at r s ~ 1 before going to zero for r s — > 0. 

The results obtained previously by Wang and Ma— for Z = 1 are qualitatively similar to our results although 
there are some quantitative differences. We recall that they approximated the screened potential by the linearized TF 
potential and adjusted a screening parameter in order to satisfy the FSR. This potential is everywhere attractive and 
misses the repulsive region shown in Fig. 1 which is associated with H~ formation at low densities. These differences 
become more apparent when the stopping powers for Z = 1 and Z = — 1 are compared. This is shown in the inset to 
Fig. 11 where the ratio of these stopping powers is given. The ratio tends to unity for r s — > since linear response 
is valid in this limit. However, the ratio also tends to unity for large r s since as we have seen, the H _ configuration 
is effectively the same as an external negative charge insofar as the scattering of continuum states is concerned (see 
Fig. 6). In contrast, for an assumed model potential, Wang and Ma find the ratio of stopping powers to be close to 2 
for all r s . This points to the danger of using model potentials which do not properly account for the true nature of the 
electronic screening. Another difficulty of the model potential approach is that differences associated with different 
screening approximations are no longer apparent. This is emphasized in Fig. 11 by the plot of the stopping power in 
the Hartree approximation. The results are qualitatively similar to those with xc, but there are important differences. 
Clearly additional input is needed in a model potential approach to capture these differences. 

Finally, we compare with the results of Krakovsky and Percus 2 ^ for Z = — 1. As mentioned in the Introduction, they 
obtain the screening charge density using the Sjolander-Stott integral equation 22 which is a different way of introducing 
exhange and correlation. This induced density is shown in a second paper 3E > and is qualitatively similar to our results 
in Fig. 5. They then use this density to construct an effective scattering potential for which the scattering phase shifts 
are determined. However they give no details about how their potential is constructed nor what interaction effects, 
such as exchange and correlation, are included. Nevertheless, the stopping powers they obtain for Z = — 1 in 2D are 
similar to our results. We find in fact that our Hartree and xc stopping powers bracket their results with deviations 
which are typically less than 20%. 

IV. HIGH DENSITY LIMIT 

We briefly commented earlier on the behaviour of the screening cloud in the high density limit (E F — > 00, r s — > 0). 
In this Section we investigate this question in more detail and in particular, establish the connection between linear 
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and nonlinear screening for a 2DEG. This problem in 3D was touched on in an early paper by Butler— which dealt 
with the possibility of positronium formation in metals. Using a square well model potential, it was shown that the 
total electron density at the origin is an analytic function of the strength of the potential. This in fact is a special 
case of a more general result obtained by Kohn and Majumdar— concerning the analyticity of physical properties of 
the entire electron gas. Butler then used his result regarding the analyticity of the density to argue that the total 
density at the position of a positron in a metal could be obtained by perturbation theory (to all orders) even if the 
screened positron potential had a bound state. This assertion was confirmed for the specific example of the square 
well model where lowest order perturbation theory already provided a good estimate of the exact total density in the 
limit of high electron gas densities. 

We shall now prove a theorem which demonstrates that some of Butler's conclusions have a much broader range 
of validity. It is unnecessary to specify either the form of the externally imposed potential nor the dimensionality of 
the system. The physical situation we consider is a noninteracting Fermi gas in the presence of an external potential 
XV (r). The effect of electron-electron interactions will be dealt with later. The potential is assumed to be smooth and 
to approach zero sufficiently rapidly as |r| — > oo. The parameter A is an arbitrary coupling constant whose physical 
value is unity. 

The problem is to determine the induced density Ara(r) which is due to the introduction of XV (r). This quantity 
is conveniently defined in terms of the Green's operator 

= ihi = T^nhxv < 12 > 

where Hq is the unperturbed hamiltonian. In a spatial representation, the Green's function is given by 



G(r,r',z) = (t\6(z)\t>) 

= ^(r)^(r') 



(13) 



z-Ei 

where i/Ji(r) is an energy eigenstate of H . From this we see that 

V \Mv)\ 2 5{E - Ei) = -~ImG(r, r,E + ie) . (14) 

i 

Thus the total density is given by 

n(r)= l^(r)| 2 = -- / *" dE ImG(r, r, E + ie) , (15) 

i(occ) 

where the upper limit of the energy integration is the Fermi energy Ep below which all states are occupied. 
The change in density from the unperturbed situation is clearly 

An(r) = n(r) — no(r) 

= -- { dEIm[G(r,r,E + ie)-G {r,r,E + ie)} , (16) 

T J-oo 

where Gq(t,t,z) is the Green's function for A = 0. It is the screening density Ara(r) that we wish to investigate in 
the high density limit. Since the potential is multiplied by the coupling constant A, the density as given by Eq. 11 tip 
depends implicitly on this parameter. The variation of the induced density with A, holding Ep fixed, is simply 

8An(r-X) 1 f E \ . .dG{E + ie) l 

Im(r| ^— -|r). (17) 



dX it 7_oo dX 

By making use of the Dyson equation 

G(z, X + SX) = G{z, X) + G(z, X)SXVG(z, X + SX), (18) 

we find 

9G d\ X) =G(z,\)VG(z,\). (19) 
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Thus, the variation of the induced density can be expressed in the form 
dAn{r;X) 1 r r^F 



OX 



= -- f dr'V(r') [ F dElm[G(r,r' , E + ie)G(r' ,r, E + ie)] . (20) 

7T J J-oo 



We note at this point that this quantity still depends implicitly on A through the full Green's function G(r,r', z). 

We now make use of the spectral representation of the Green's function in Eq. (|13|l . The product tpiir)^* (r 1 ) can 
be treated as real since a complex wave function at the energy Ei must have a time-reversed pair which is its complex 
conjugate. Thus, we have 



5An(r; A) 



dv' V(r') V MrMir'^ir'Wjir) 



d\ J v ' ^ 

y 

i rEp 

x P- / dE [5{E - Bk) - 5(E - E } )] . (21) 

This expression is formally exact and is the crucial step in the argument needed to arrive at the final conclusion. It is 
clear that the only terms which give a nonvanishing contribution are those with Ei < Ep and Ej > Ep, or vice- versa. 
It is these terms which give a finite result in the Ep — > oo limit. If this limit is taken first for a given (i,j), the result 
is zero. In other words, it is important to sum over all states before the Ep — > oo limit is taken. We also note that 
the i = j terms in Eq. (|21|l do not present any difficulty despite the appearance of the energy denominator. From 
Eq. H2U|) . these terms contribute 



^|^(r)| 2 |^(r')| 2 Imhm J^dE ^-1 

dE 

-oo 



-Ei) 




2 



dE E + ie- E, 

Tr^lM^mirTSiEp-E,), (22) 



which is well-defined. 

The energy integral in Eq. H21[) can be written alternatively as 



/Ep r oo 

dE[S(E - Ei) - S(E - Ej)} = - / dE [5(E - - S(E - Ej)] . (23) 
-oo J Ep 



Using this in Eq. J2IJ, we find that Eq. Q20(l is equivalent to 
<9An(r;A) _ 1 



OX 



= - f dr'V(r') f dETm[G(r,r',E + ie)G(T',T,E + ie)] . (24) 

n J JEp 



This is a more useful form for taking the Ep — > oo limit since now the energy argument of the Green's function is large. 
If Ep ^> | V(r) | everywhere (we assume here a bounded potential—), the effect of the potential on the propagation 
of the electron can be neglected and the Green's function G can be replaced by the free-particle propagator Go- We 
then have 



dX 



dAn{r;X) _ 1 f dj jy^ Hm f dEIm[G(r,r' , E + ie)G(r' ,r, E + ie)] 

7T J E f ~,ooJ Ef 

~ - [ dr'V(r') lim / dEIm [G (r, r', E + ie)G (r', r, E + ie)] 

7T J E f-'°°Jep 

= -- [ dr'V(r') lim ( dE Im [G (r, r', E + ie)G (r', r, E + ie)] , (25) 

7T J B F-»°°7-oo 



where we have simply reversed the earlier argument to change the limits of integration in the last line. Several 
comments are in order: 

1. The replacement of G by Go is only valid if the energy argument is large; these Green's functions are certainly 
different at low energies. In fact, the perturbed and unperturbed systems have completely different state spectra 
and H may contain bound states which do not appear in Hq. Nevertheless, since Eq. (|20|l can be written in the 
equivalent form (|24|l . the replacement of G by Go in Eq. (|25|l is justified in the Ep — > oo limit. 
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2. The final line in Eq. (|25() is identical to the result obtained by linear response theory. Since it is independent of 
A, it can be integrated from A = to A = 1 to give 



An(r) = -- I dv'V{r') lim / dE Im [G (r, r', E + ie)G (r' , r, E + ie)] 



(26) 



3. The final result in Eq. (|26|l does not depend on the assumed dimensionality and is valid even when the potential 
V(r) supports bound states. These bound states of course do not appear in the spectrum of Go but are 
nevertheless included in An(r). The result given in Eq. (|26|l will be referred to as the high density screening 
theorem. 

As an application of Eq. (12(:i[) we consider an ideal uniform gas in D dimensions. The use of the spectral represen- 
tation for the Green's function gives in this case 



- dE lm.[G (r,r',E + ie)G (r',r,E + ie)] 

K J-oo 

kk' k k 

i_^ e ^(r-r') Xo(q) 



L D 



(27) 



where the factor of 2 is for spin, L is the size of the sample in a given direction, f(E) is the zero-temperature Fermi 
function and Xo(*l) is the static free particle density response function 



/(£ k+q ) - f(E k ) 



Substitution of Eq. J23 into Eq. igEJl gives 

An(r) = - lim -L £ Xo(q)V(q)e* 



(28) 



(29) 



This is a rigorous statement that linear response gives the correct induced density in the high density limit in any 
dimension. 

The static response functions in the various dimensions arc given by the expressions 



2k F + q 



2k F — q 



xl D (q) = -e(2k F -q) + - 1 - v/1 - {2k F /qf 6{q-2k F ) 

IT IT 

k F l fcf 



2tt 2 



1 
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1 



g 
44 



In 



2fci? — g 



2k F + q 



In each case, the q — > limit of xo(<?) is the density of states g(E F ) at the Fermi level 

2 



Xo» 



1 

7T 
A'/: 



(30) 



(31) 



If the potential is smooth we can suppose that V(q) is effectively zero for q > q max with q max <C k F . In this case, 
Eq. reduces to 



An(r) = - lim xo(0)V(r), 



(32) 
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which is recognized as the linearized TF approximation for the density. 

It should be noted that 2D is special for a number of reasons. First, since Ep is proportional to the density n, the 
full TF approximation is itself linear. Second, the free-particle response function has the constant value Xo(q) = 1/tt 
for all wavevectors up to 2kp. Thus, the replacement of Xo(q) by Xo(0) is in some sense a less restrictive approximation 
than it is in other dimensions. In particular, if the only finite Fourier components of V(r) occur for q < 2kp (which 
is always the case when kp — * oo), we have 

Ara(r) = — V(r) . (33) 

The analogous result in 3D is An(r) = while the ID result does not have a finite limit. Thus the induced density 
only has a nontrivial limiting value in 2D. 

The result for the 2D screening charge in Eq. i|33|) is all the more remarkable if we note that the potential V(r) can 
in fact support bound states, as it must in 2D if the potential is purely attractive. It is worth verifying this explicitly 
by means of a numerical example. We consider the following model potential 

V{r) = Sv s^(^), r<r , (34) 
{ , r > r , 

which has no physical basis but is chosen simply to illustrate the point. It is a repulsive double barrier potential for 
Vo > and an attractive double well potential for Vb < 0. With this potential we calculate the induced density by 
solving the Schrodinger equation for all states below Ep. In Fig. 12 we show results for a repulsive potential with 
Vo = 0.125 H and r — 5 a.u. as a function of gas density. The screening charge is due to continuum states in this 
case and it can be seen that by r s = 0.5, the induced density is very close to the prediction of Eq. (|3"3T) . Similarly, if 
an attractive potential is placed in a gas with r s — 0.5 (Ep — 4.0 H), we obtain the bound and continuum charge 
densities shown in Fig. 13. Here, the magnitude of \Vq\ was adjusted to vary the number of bound states. Two cases 
were considered: (a) Vq = —0.125 H for which a single m — bound state is found and (b) Vq = —0.25 H, for which 
two m — bound states exist. The bound and continuum densities are quite different for these two cases, but the 
total density agrees very well with the result in Eq. 1)33(1 , even though the high density limit has only marginally been 
reached. These examples confirm very nicely the behaviour stipulated by the high density screening theorem. 

We can now apply the theorem to explain some of the results found in the Z — 1 calculations of Sec. III. Of course 
interactions were included there, but the Kohn-Sham system of equations correspond to a system of noninteracting 
electrons moving in the total effective potential in Eq. @. In Fourier space, we have 

Ai>eff(q) = v ext (q) + v(q)An(q) + v' xc (n Q )An(q) . (35) 

We have here linearized the xc potential since, in the high density limit of interest, the induced density is small 
compared to the background density. Identifying Av c s(q) as given by Eq. I|35(l with V(q) in Eq. (|29|l . we have 

An(g) = - X o(q)Av eS (q) (36) 

which implies 

An / a = Xo(gKxt(<?) , 37 x 

1 + (v(q) + V xc (n ))xo{q) ' 

This is simply the induced density obtained within linear response theory, treating Hartree and xc interactions at 
the mean-field level. An example of the linear response result was already given in Fig. 9 for the case of a Z = 1 
charge outside the plane of the 2DEG. The close agreement of this density with the nonlinear screening density is 
confirmation of the high density screening theorem. 

For high densities, v xc is in fact dominated by the exchange potential v x = —C x /r s and we have v' xc (rio) ~ 
— (irC x /2)r s . Thus in the limit r s — > 0, the xc term can be dropped and the induced density is simply given by the 
result 

i + «(?)xo(0) 

The corresponding effective potential is 



(39) 
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The real space quantities are obtained from these expressions by an inverse Fourier transform. 

We now consider the particular example of a point charge Z in the plane of the gas. With v(q) — lir/q in 2D, 
Eq. (J3HI) becomes 



A «cff(<z) = -^, (40) 



and in real space we obtain the TF potential 8 



"tfM = -Z I dq—^-J (qr) 



q + 2 



= -z[ ^J^W < 41) 

It is a monotonic potential which behaves as —Z/r for r — > and falls off as — Z/4r 3 for r — > oo. Being negative 
definite, it will support at least one bound state for any value of Z. For Z = 1 the bound state energy is found to be 
E = —0.2862 H which agrees with the value found by Stern and Howard 8 -. This is the limiting value of the bound 
state eigenvalue shown in Fig. 3. 

We thus come to the conclusion that the r s — > limit of the Kohn-Sham effective potential is the TF potential in 
Eq. glj). To illustrate this we show in Fig. 14 the induced density for Z = 1 and a very high density corresponding to 
r s = 0.2. Also shown are the bound and continuum components which add to the total. The latter is also displayed 
as rAn(r) in Fig. 15 together with the TF density as obtained from Eq. I|38|) . It can be seen that the DFT density 
simply oscillates about this limiting value. With decreasing r a , the amplitude of these oscillations decrease and the 
DFT density uniformly approaches the TF density. 

For completeness, we note that the result corresponding to Eq. ijlUI) in 3D is 

Av eS ( q ) = - -2— 5 , (42) 



with qxF = y/^kp/w the Thomas- Fermi wavevector. The real-space potential in this case is 

Ze~ qTFT 

v eS (r) = . (43) 

r 

With increasing hp, the range of the potential decreases. Thus at some point the potential will no longer support 
a bound state, and all states in the presence of the external charge will be continuum scattering states. As alluded 
to earlier, the screening of the external potential becomes more effective with increasing density. This is in marked 
contrast to the situation in 2D where the screened potential has a density independent limit which continues to 
support bound states in the high density limit. 

In summary, if a potential V(r) is introduced into a 2D gas and the induced density is constructed by summing 
up over all eigenstates of the hamiltonian H (including all possible bound states), the induced density will be given 
exactly by Eq. 1)33(1 in the high density limit. Remarkably, this density is the same as what one would obtain by 
treating the potential as a weak perturbation and applying perturbation theory. 



V. CONCLUSIONS 



In this paper we have provided a detailed discussion of the nonlinear screening of an external point charge by a 
two-dimensional electron gas. Nonlinear effects are most evident when the charge lies in the plane of the 2DEG 
but rapidly decrease in importance as the charge is moved out of the plane. As an application of the theory we 
have considered the problem of 2D stopping power and have compared our results with earlier work. Quantitative 
differences point to the importance of performing self-consistent nonlinear screening calculations. 

We also derive what we refer to as the high density screening theorem. In 2D, the theorem has the interesting 
consequence that the screened potential has a density independent limiting form in the high density limit and that the 
screening density is directly proportional to this potential. Explicit calculations for model potentials confirm these 
conclusions and demonstrate the surprising fact that linear response theory is valid even when the potential acting on 
the gas supports bound states. These nonintuitive results highlight some of the peculiarities of electronic screening 
in two dimensions. 
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APPENDIX A: EVALUATION OF HARTREE POTENTIAL 



In this Appendix we summarize the methods we have used to determine the electrostatic potential of the screened 
impurity charge. Due to the long-range nature of the Coulomb interaction, this is the one nontrivial part of the 
calculation of the total self-consistent potential. 

A charge Z at the position R = dz above the plane of the 2DEG [z = 0) is represented by the external charge 
density 

n oxt (r) = Z6(r-K) (Al) 

which gives rise to the electrostatic potential 

^xt(p) = {p2 + d2)1/2 (A2) 

in the plane of the gas. Here, p is the radial distance from the origin. This potential induces the electronic screening 
density An(r) = An(p)S(z) which in turn gives rise to an electrostatic potential that is determined by Poisson's 
equation 

V 2 A0(r) = 47rAn(r) . (A3) 

Although this equation determines A<j) in three-dimensional space, only the behaviour within the plane of the gas is 
required in Eq. 

To solve Poisson's equation, we make use of the 2D Fourier transform pair: 



A4>(q,z) = I Spe-^PA^z) 

(27T) 2 



A <M^) = / 77dvl e iq ' P A#7^) • (A4) 



Taking the Fourier transform of Eq. (|A3|) , we obtain 



dz 2 



g 2 A0 = 4irAfi5{z) , (A5) 



which has the solution 



with 



A4>(q, z) = A4>(q, z = 0)e- q ^ (A6) 



A4>(q,z = 0) = -— An. (A7) 

q 

The screening potential in the plane of the gas is then obtained by taking the inverse 2D Fourier transform of this 
equation. 

Two problems arise in the numerical implementation of this procedure. First, Eq. (|A7|I has a q = singularity 
associated with the long-range nature of the Coulomb potential, and second, a straightforward iteration of the Kohn- 
Sham equations does not converge. To deal with the first problem, we introduce an auxiliary charge density n aux (p) 
which gives rise to the auxiliary electrostatic potential 

aux (p) = {p2 + Z dl)1/2 , (A8) 
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where do is some fixed, nonzero parameter. Evidently, this potential can be considered as being due to a point charge 
at R = doz. Alternatively, its 2D Fourier transform 

kux(q) = ^-Ze-" d0 (A9) 

implies that the auxiliary charge density in the z = plane has the 2D Fourier transform 

n aux (q) = Ze-^ . (A10) 

The corresponding real space density is therefore 

n * M = ^ { p i + dl)m - ( An ) 

This density has a total charge of Z and falls off as p~ 3 as p — > oo. 
We now write the external potential as 

0ext(p) = 0aux(p) + A<j> cxt (p) (A12) 

where 

1 1 



A0 cx t(p) = Z 



(p 2 + d 2 )V2 (p2 + d 2y/2 



(A13) 



The latter is treated explicitly in real space, while the effect of (/> aU x is accounted for by means of the following Poisson 
equation: 



d 2 Act) 



dz 2 

This equation gives 



q 2 A<j> = -4re-(n aux - An)S(z) , (A14) 



2ir 

A(f>(q, z = 0) = — (n awc - An) (A15) 

q 

which no longer has ag = singularity since the combination fi aux — An is charge neutral. Thus, its inverse Fourier 
transform can be calculated readily using a fast Fourier transform (FFT) technique, and the addition of the result to 
Eq. i|A13(l yields the total electrostatic potential. 

The second problem concerns the convergence of the self-consistent iterative procedure. This is dealt with by 
rewriting Eq. I|A14(I in the form 

- (q 2 + * 2 )A4> - -47r(n aux - Ah)S(z) ~ k 2 A4> . (A16) 

This trivial change obviously does not change the solution of the equation but we see that the left hand side corresponds 
to a screened Coulomb interaction with screening parameter k. With this equation in mind, we adopt the following 
iterative procedure^: 

A~^ +1 \q,z = 0) = ^(n aux - An«) + (l - ±\ A~^{q,z = 0) , (A17) 

where Q 2 = q 2 + n 2 . It is evident that when A0^ +1 ^ = A^W at convergence, this procedure yields 

A4>(q, z = 0) = y (n aux - An) , (A18) 

which is the correct potential for the charge density n aux — Are. The use of a screened- Coulomb interaction has the 
important effect of stabilizing the iterative procedure required to achieve selfconsistency. 
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FIG. 1: Self-consistent effective potentials Av e s(r) in Hartree units as a function of position. Each solid curve corresponds to 
a different r s value, ranging from 1 for the lowest curve to 10 for the highest curve, in steps of 1. The dashed curve is the 
Thomas-Fermi screened potential. 
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FIG. 2: Screening density as a function of position for two r s values. The dashed curve is the bound state contribution, the 
chain curve is the continuum state contribution and the solid curve is the total. 
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FIG. 3: Bound state eigenvalue in Hartree units vs. r 3 . The solid curve is for the full nonlinear DFT calculation including 
exchange-correlation, while the dashed curve is the corresponding result in the Hartree approximation. 
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FIG. 4: Self-consistent effective potentials with (solid curve) and without (dashed curve) exchange-correlation. 
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FIG. 6: Normalized density as a function of position for r s = 10. The solid curve is for Z — 1 while the dashed curve is for 
Z = -l. 
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FIG. 7: Comparison of nonlinearly screened and linearly screened potentials for r s = 4 and Z — —1. The solid curve is 
the nonlinearly screened potential including exchange-correlation, u ex t + Avh + Aw^ c ; the short-dashed curve is the Hartree 
potential, v cxt + Avh, calculated with the same screening density as used for the solid curve ; the chain curve is the linearly 
screened potential within the RPA (i.e., the Hartree approximation); and the long-dashed curve is the Hartree potential, 
Uoxt + Avh, obtained with a linear screening density that includes local-field effects (i.e., exchange-correlation). 
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FIG. 8: Screening charge density as a function of position for Z = 1 and r 3 — 2 for different distances d of the charge from the 
plane of the electrons: d = 0.5ao (solid curve), d — 1.0a o (chain curve), d — 2.0a ( dashed curve). 
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FIG. 9: Screening charge density as a function of position for Z = 1, r a = 2 and d — 2.0cio. The chain curve is the bound state 
contribution, the long-dashed curve is the continuum contribution and the solid curve is the total. The short-dashed curve 
almost coincident with the solid curve is the linear response result including exchange-correlation. 
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FIG. 10: Normalized stopping power as a function of the projectile charge Z, for r s — 2 and d = 0. The solid curve is 
obtained from the nonlinear screening result including exchange-correlation; the dashed curve shows the corresponding result 
in the Hartree approximation (no exchange-correlation); the straight line (chain curve) is the quadratic response result— 
{S 1 +S 2 Z)/v. 
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FIG. 11: Stopping power as a function of density parameter r s for Z = — 1. The solid curve shows the full nonlinear screening 
result, including exchange-correlation, while the dashed curve gives the corresponding result in the Hartree approximation. 
The inset shows the stopping power ratio S(Z = 1)/S(Z = —1) as a function of r s . 
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FIG. 12: Screening charge density as a function of position for the model potential in Eq. 1341 . The potential parameters are 
V = 0.125 H and r = 5a . 
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FIG. 13: Screening charge density for r s = 0.5 for the model potenial in Eq. I|34|l : (a) Vo ~ — 0.125 H, ro — 5ao; (b) 
Vo = —0.25 H, ro = 5 ao. The dashed curve is the bound state contribution, the chain curve is the continuum state contribution 
and the solid curve is the total. 
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FIG. 15: As in Fig. 14, but plotting the screening charge density as rAn(r) vs. r. The dashed curve is the limiting (r s — > 0) 
Thomas- Fermi density. 



